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Abstract 

A consistent and efficient set of boundary 
conditions are developed for the multi-sweep 
space marching pressure-elliptic Reduced Navier- 
Stokes (RNS) scheme as applied for three- 
dimensional internal viscous flow problems. 
No-slip boundary conditions are directly imposed 
on the solid walls. There is no iteration proce- 
dure required in the cross plane to ensure mass 
conservation across each marching plane. The 
finite difference equations forming the coeffi- 
ceint matrix are ordered such that the surface 
normal velocity is specified on all the solid 
walls; unlike external flows, a pressure boundary 
condition in the cross plane is not required. 
Since continuity is directly satisfied at all 
points in the flow domain, the first order momen- 
tum equations can be solved directly for the 
pressure without the need for a Poisson pressure 
correction equation. The procedure developed 
herein can also be applied with periodic bound- 
ary conditions. The analysis is given for gen- 
eral compressible flows. Incompressible flow 
solutions are obtained, for straight and curved 
ducts of square cross section, to validate the 
procedure. The solutions of these test cases 
are used to demonstrate the applicability of the 
RNS scheme, with the improved boundary condi- 
tions, for internal flows with strong interac- 
tion as would be encountered in ducts and 
turbomachinery geometries. 

Introduction 

The flow through advanced highly loaded 
turbomachinery blade rows is characterized by 
extensive regions of strong three-dimensional 
viscous-inviscid interaction. To simulate this 
important pressure interaction, numerical meth- 
ods that can couple the viscous and inviscid 
regions must be employed. In addition, effi- 
ciency and accuracy of the numerical algorithm 
become important considerations if these methods 
are to be useful in the aerodynamic design 
process. 

A number of different approaches have been 
considered for the simulation of strong viscous- 
inviscid interaction in internal flows. Conven- 
tional full Navier-Stokes (N-S) methods, which 
solve the full N-S equations throughout the flow 
field, have been successfully used to analyze 
three-dimensional interacting flows in turboma- 
chinery blade passages. 1 * 2 However, these meth- 
ods do not exploit the asymptotic behavior of 
the equations at the large Reynolds numbers typi- 


cally encountered in turbomachinery flows. Con- 
sequently, these methods require large computer 
storage and run times. A recently developed 
method in this category 1 requires 18 hr of CPU 
time on a VAX 11/780 computer for the prediction 
of end-wall flow in a cascade on a relatively 
coarse grid of 53 by 31 by 10 nodes. 

Another approach is interacting boundary 
layer theory. This has been used by a number of 
researchers 2-7 for two-dimensional applications, 
where the interaction of the inviscid flow on the 
boundary layer is coupled through the injection 
and surface boundary conditions for the inviscid 
and boundary layer analysis respectively. These 
methods are potentially very efficient; however, 
the evaluation of the injection condition or 
inviscid displacement body (due to the viscous 
effects), which alters the inviscid flow, can 
become rather involved for complex three- 
dimensional flows. Approximate methods normally 
used to evaluate this effect, such as linearized 
small disturbance theory, can result in consider- 
able error. In addition, the approximation of 
zero normal pressure gradient through the bound- 
ary layer might not be appropriate for turbulent 
flows with strong pressure i nteraction . 1 6 

Methods that use space marching with an 
approximate form of the steady N-S equations 
(single-pass and multipass marching methods) 
have been considered for a number of years to 
predict flows through curved ducts and turboma- 
chinery blade cascades . 8-1 1 Single pass march- 
ing can be used for configurations where the 
flow is of initial value character, but multiple 
pass procedures are required for elliptic flows. 
When applied to elliptic flows, these formula- 
tions have generally introduced a Poisson equa- 
tion for pressure to correct an initially 
assumed pressure field. This is required in 
lieu of the continuity equation, which is not 
satisfied explicitly, in order to ensure global 
mass conservation. These methods are called par- 
tially parabolic or semi-elliptic methods to 
distinguish them from the full N-S schemes. 
Although these methods result in less computing 
time than full N-S methods, the solution of the 
Poisson equation still requires large computer 
run times. In addition, due to the uncoupled 
nature of the pressure correction, which is a 
necessity of the formulation, extremely large 
under-relaxation is required in high subsonic and 
transonic flow regions. This slows down the con- 
vergence of the iteration procedure and thereby 
further increases the run time. 
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A method that combines the asymptotic treat- 
ment of interacting boundary layer theory and the 
accurate interaction simulation of the full 
Navi er-Stokes methods is the Reduced Navier- 
Stokes (RNS) formulation. This scheme was origi- 
nally developed for external flows 12 * 1 ** anc j later 
formulated for internal flows. 19 The solution 
procedure and the boundary conditions have been 
modified in this study to make the scheme more 
efficient for both two- and three-dimensional 
flows with strong interaction. As described in 
earlier references , 14-19 the system of equations 
resulting in the RNS formulation is similar to 
that of the partially parabolic scheme in that 
streamwise diffusion effects are neglected. How- 
ever, the elliptic effect or upstream influence 
in strongly interacting flows is simulated by a 
characteristic treatment of the streamwise pres- 
sure gradient. The solution procedure is there- 
fore very much different, and more direct, as 
compared to that of partially parabolic schemes. 
The equations are solved by a relaxation proce- 
dure with full coupling between pressure and 
velocities and without the need for a Poisson 
equation for the pressure correction. Detailed 
analysis of the RNS scheme and solutions for lam- 
inar, turbulent, subsonic, transonic, and super- 
sonic flow regimes for a variety of external 
flow configurations are given in Refs. 12 to 18. 
Application of the scheme for internal flow and 
some preliminary results for two- and three- 
dimensional internal flows were presented in 
Ref. 19. As pointed out in these references, the 
procedure is applicable to both inviscid and vis- 
cous flow and can be classified somewhere between 
interacting boundary layer theory and full 
Navier-Stokes solvers. 

A detailed description and analysis of the 
RNS scheme as applicable to two-dimensional 
external flows are given in Refs. 13 to 18. 
Details of various stages of the evolution of 
the scheme leading to its present form are also 
given in Refs. 13, 15 to 17. For the sake of 
completeness, some of the analysis is repeated 
here . 

The RNS equations were first considered as 
single sweep or PNS (Parabolized Navier-Stokes) 
marching procedures for supersonic flows. The 
first application was for hypersonic flow 20 
where the contribution of the streamwise pres- 
sure gradient in the corresponding momentum 
equation is negligible and can therefore be 
neglected. The equations are mathematically par- 
abolic, upstream influence is negligible, and an 
exact solution is obtained in a single marching 
sweep. For lower supersonic mach numbers, where 
the influence of the streamwise pressure gradient 
is not negligible, an elliptic effect associated 
with pressure interaction through the subsonic 
portion of the boundary layer introduces upstream 
influence. 21 ’ 22 A single sweep methodology then 
leads to an ill-posed initial value problem and 
gives rise to exponentially growing departure 
solutions for a marching step size, A£ , less 
than (A£>min where (A ^min 1s proportional 
to the extent of the subsonic portion of the 
flow in the normal (cross stream) direction. 14 ’ 23 
To surpress this so-called departure effect that 
reflects the boundary value character of the 
problem, researchers have used a variety of 
approximation techni ques 2 ^ -24 to simulate the 


elliptic effect of the streamwise pressure gradi- 
ent term. 

In the present RNS procedure, the streamwise 
pressure gradient term is split according to 

its characteristic behavior so that 

= u(P^h + (1 - u)(P$r)e 

This follows the eigenvalue analysis of Vigneron 
et. al, 25 where 0 < w(M) < u max is a function of 
local Mach number M and 

“max - {yM 2 /(U( Y - 1 >m 2>, l} min 

As mentioned in Refs. 18 and 19, the por- 
tion o)( p£)h» which is "backward” differenced 
during discretization, represents the "hyper- 
bolic" or marching part of Pr and the term 
(1 - w)(P^) e represents the "elliptic" or 
relaxation contribution that is "forward" differ- 
enced. Note that for i mcompress i bl e flow, 
since to(0) • 0, the entire P^ contribution 
iselliptic. "Forward" differencing of 
(1 - u)(Pc) e introduces upstream influence in 
the computational domain. This removes the 
i 1 1 -posedness found in the single sweep initial 
value formulation. Due to the forward differenc- 
ing, the solution procedure requires multiple 
sweep marching or relaxation. The above treat- 
ment of the streamwise pressure gradient, with 
multiple sweep relaxation, leads to consistent 
(arbitrary A£) and departure free ( A£ ■* 0) solu- 
tions for the entire range of incompressible to 
supersonic Mach numbers. Significantly, only the 
pressure (and possibly the axial velocity in the 
limited regions of reversed flow only) need be 
stored. This results in, among other advantages, 
a significant reduction in storage requirement 
over conventional N-S methods. 

In the present study, a new consistent solu- 
tion procedure is formulated, using the RNS 
Scheme, for three- and two-dimensional flow prob- 
lems. The treatment of boundary conditions has 
been significantly modified, compared to the ear- 
lier procedure, 19 to make the solution procedure 
more efficient and accurate. The application of 
zero injection or solid wall boundary conditions 
in the cross-plane is more direct in this study 
than in that of Ref. 19. The procedure is deve- 
loped for arbitrary compressible flow, but only 
incompressible solutions are obtained for deve- 
loping flow in three-dimensional straight and 
curved ducts of square cross section. These 
solutions are compared with available experimen- 
tal data and computed results. 


Governing Equations 


The governing equations are written in a gen- 
eral curvilinear coordinate system (£, n, and £) 
in terms of the primitive variables (u, v, w, p). 
The momentum equations are then rearranged to 
reflect the momentum balance in the directions of 
the contravariant velocity components (u, V, 
and W). This requires the appropriate combina- 
tion of the Cartesian component momentum equa- 
tions after transformation into the £, n, and 
£ coordinate system. For example, the momentum 
equation in the U direction is written as ? x 
(x-momentum) + £y(y-momentum) + £ z (z-momentum) . 2 *> 
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The final equations are given in the following 
matrix form. 


Continuity and Momentum: 


3jrE + 9pF + 9^G = 3^R + + 9^T + K 

where 


E = 3 


'pCU ] 
P^uu 
pCUV 
-p;uwJ 


[pCV 1 
pCVU 
pCW 

LpCVwJ 

fpCW ] 
pCWU 
pCWV 

Lpcwwj 


E xV 


Vx 


Vx 


+ Vy 5 + c z T z e 
c + Vy C + Vz C 
; + Vy ? + Vz C 


Vx 


Vx 

Vx' 


+ Vy + Vz 
+ Vy" + Vz n 
+ V y n + Vz n j 


W 


E xV 


Vx 

Vx 


+ Vy C + Vz ? 


+ n y T y + Vz 
+ V/ + Vz 


K = 3 


p^G 


« 


p;g 


u 


pCG 


« 


- - g 5n P n - g 5 S|- 

- - g 5n P n - g n? P c 

- g 5? p ? - g^ n P n - g c S ? 


The terms x S, x 5, etc. appearing in the 
x y 

column vectors R, S, T, and K are explained 
in detail in Appendix A. The velocities U, V, 
and W are the contravariant components; all 
the shear stresses, as shown as Appendix A, can 
be expressed in terms of these components. 

Since one of the coordinates (f-) represents the 
marching direction, it has been found that the 
equations in this form enhances the stability of 
the numerical scheme. In addition, the system of 
equations in this form can be easily verified 
for an orthogonal coordinate system. For the 
sample problems considered in this study, an 
orthogonal (curvel i near) coordinate system is 


specified. For the present formulation, adia- 
batic (wall) conditions are assumed and with a 
Prandtl number of unity, a simplified energy 
equation results; i.e., total enthalpy is con- 
stant. The same algorithm can also be used for 
nonadiabatic wall conditions and Prandtl numbers 
different from unity. The energy equation, writ- 
ten in terms of stagnation enthalpy is only 
weakly coupled with the remainder of the equa- 
tions for low speed and even moderate supersonic 
flows. Therefore, the energy equation can be 
solved in an uncoupled manner to update the stag- 
nation enthalpy at each marching location. To 
close the system, an equation of state p = pRT 
and a relation for viscosity, y = y(T), are 
required for compressible flows. 

In all of the momentum equations, the diffu- 
sion terms in the streamwise direction are 
neglected according to the RNS approximation. 
These terms are negligible for the coordinate 
system specified herein so that the RNS system 
closely approximates the full N-S equations. 

The flow Reynolds number is based on the inlet 
uniform velocity and the inlet hydraulic radius 
of the duct. The pressure is nondimensional i zed 
with the inlet dynamic head. 

Pi scretization 


The discretization of the governing equa- 
tions is illustrated in Figs. 1 and 2. In the 
marching (£) direction the equations are discre- 
tized at (i), the velocity node points. In the 
cross-plane, the continuity, streamwise momentum 
and the two normal momentum equations are discre- 
tized at O,(3)>0> and (p, respectively. All 
of the connective terms in the marching direc- 
tion are upwind differenced. Both first-order 
(two point) and second-order (three point) accu- 
rate upwind differencing schemes are considered. 
As discussed earlier, the streamwise pressure 
gradient, P$r is differenced as: 


P 5 " (1 - “1+1 5 ( 


,n-l 


ni 


HI - Pi 


5ui - E 1 


+ (/ V f 


r n n 
P 1 - P 1-1 


Vni j 


+ <Ap e } bl 


where the subscript i is a modified index for 
pressure in the ^-direction (see Fig. 1) and 
n is the current marching (global) sweep. The 
terms (Ap^ and (Ap^ are additional correc- 
tion terms to produce second-order accuracy 
in the forward and backward directions respec- 
tively. These terms are given by: 

<APE> f = a f (l ♦ a f )<C 1+1 - Cj> 


X 


^( 1 +Of) 




■/] 


*3 



and 


Boundary Conditions 


(Ap?>b = v 1 + V <c i - «i-i> 


x [o b p" - ( 1 + ° b > l>”_ + pj 


Of = <51 + 1 - 5i)/C?i - Ci-1) 

*b - <51-1 - * Si-1) 

For first order accuracy the terms (Ap$r) f and 
( Ap$r ) ^ are neglected. 

The discretization for p requires that 
the unknown pressure p. at the marching loca- 
tion i, is staggered at a distance (l-w)A£ 
upstream of velocity Uj , , and Wq . The 

pressure at the grid point i is given by: 

Pi = “i + l p i + (1 - “l + l^l+l 

The discretization of all derivatives in the n 
and t; cross flow directions is second-order 
accurate except for one of the viscous terms in 
each of the normal momentum equations. These 
terms are the diffusion terms in the same direc- 
tion as the momentum direction, i.e., the second 
derivative with respect to n in the n-momentum 
equation and the corresponding term in £ momen- 
tum equation. Since these viscous terms are 
negligibly small, i.e., of order of the neglected 
streamwise diffusion terms, for large Reynolds 
numbers, the discretization in the n and £ 
directions remains essentially second-order accu- 
rate. The nonlinear convective terms in the 
finite difference equations are qual si 1 i eari zed , 
with respect to the previous marching location, 
using a Newton linearization. The linearized 
equation result in a coupled system for the vec- 
tor [U, V, w, p] T . The system of equation can 
be represented by the matrix equation 

[A] {$} = {q} 

where {$} is the solution vector [U,V,W,p] T , 

{q} is the known right hand side of the equa- 
tion, and [A] is the nine diagonal coefficient 
matrix shown in Fig. 3. The associated computa- 
tional molecule is shown in Fig. 4. Each of the 
elements of the coefficient matrix [A] is a 4 
by 4 matrix corresponding to the column size [4] 
of the vector {<(>}. The discrete system of equa- 
tions is solved for [U,V,W,p] T i i ^ at each 
marching location i with the modifed strong 
implicit procedure (MSIP) of Ref. 27. This 
method was originally developed for the scalar 
system describing the two-dimensional heat con- 
duction equation. 27 It was modified by the 
present investigators for application to a vec- 
tor system of equations in a previous study . 1 9 
The details of the procedure are given in 
Refs. 19 and 27. As mentioned in Ref. 19 for a 
two-dimensional flow the cross-plane reduces to 
a line and the system of equations reduces to a 
block tridiagonal system for [U,V,p]^j which 
is solved by standard LU decomposition. 


The method of application of boundary condi- 
tions often dictates the efficiency of a solu- 
tion algorithm. A major change in the RNS/MSIP 
solution procedure has been implemented in the 
present study by a modification of the boundary 
conditions as applied in the previous study. 19 

Inflow and outflow boundary conditions are 
straight forward. Since streamwise <£) diffu- 
sion terms are neglected and streamwise convec- 
tion terms are upwind differenced, the velocities 
have an initial value character (except in 
regions of reversed flow). Therefore the veloci- 
ties must be specified only at the inflow boun- 
dary assuming that flow reversal does not occur 
at the outflow boundary. Due to the splitting 
of the streamwise pressure gradient into a for- 
ward differenced (relaxation) and backward diff- 
erenced (marching) elements, a pressure condi- 
tion is required at both the inflow and outflow 
boundaries for mach numbers 0 < M < 1. For 
mach numbers M > 1 , to becomes unity and there- 
fore a pressure condition is not required at the 
outflow boundary i.e. full marching. For incom- 
pressible flow (« * 0) a pressure condition is 
not required at the inflow boundary (full relaxa- 
tion). For 0 < M < 1 , the staggered pressure 
(see Fig. 1), leads to a partially prescribed 
pressure condition at both the inflow and out- 
flow computational planes, since at the node 
poi nt 1 , 

p, = <o 1 + 1 p" + (i - « 1 + 1 )pJ;] 

In the cross plane <n - t plane), it is 
important to apply the boundary conditions in a 
consistent manner if the system of discretized 
equations is to be solved efficiently. To illus- 
trate the present procedure, let us first con- 
sider the two-dimensional flow problem. The 
cross-plane then reduces to a line, e.g., 

0 < n < 1* The discretization location for the 
continuity, ^-momentum and n-momentum, as asso- 
ciated with the node j, are denoted in Fig. 5 
by Q,0, and (rj) respectively. 

As seen in Fig. 5, the number of discrete 
^-momentum equations is j m ax - 2 ancl the n umber 
of descrete n-momentum and continuity equations 
are each j max - 1. Therefore the total number 
of unknowns is 3j max (U, V, and p for each node) 
and the total number of discrete equation is 
( 3 j max ~ 4) . For solid walls (j = 1 and 
J = jmax)* each wal1 has two Physical boundary 
conditions, i.e., U * 0 and V « 0. Therefore, 
the system is closed in so far as total number 
of equations plus the boundary conditions com- 
pared with the number of unknowns is concerned. 
However, since there are only imx - 1 discrete 
n-momentum equations, the numerical solution 
procedure apparently requires a condition for 
pressure at one of the boundaries (j * 1 or 
j * Jmax)- This would render one of the boun- 
dary conditions on V redundant. In the previ- 
ous study, 19 the zero normal velocity (V = 0) 
was indirectly imposed at the outer boundary 
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(j 55 jmax) through an artificial pressure boun- 
dary condition. An iterative process on this 
pressure boundary condition, imposed at 
j = j ma x was required in order to ensure global 
mass conservation or that the velocity V = 0 

j ■ jmax- In the present study, the zero 
injection conditions are directly imposed at 
both the boundaries without any need for the 
iterative artificial pressure boundary condi- 
tion. This is acomplished by slightly changing 
the structure of the block tridiagonal matrix 
near the outer boundary. The block tridiagonal 
system at the marching location (i), for two- 
dimensional flow, is shown in Fig. 6. For an 
interior point, 2 < j < j max _2, the three equa- 
tions for the unknown U, V, and p are 
grouped as: 

'continuity at j - 1/2' 

5-momentum at j 
.n-momentum at j + 1 / 2. 

For the point next to upper boundary (the lower 
boundary can be chosen instead of upper boundary) 
we modify the grouping as follows. 

‘continuity at (j - 1 - 1/2)' 
5-momentum at (j"~ - 1) 

continuity at (j™** - 1 + 1/2) 

For wall boundaries (i.e., j = 1 and j = j max ) 
the structure is given as follows. 

For j = 3max 

'U - 0 

V = 0 (or specified) 

n-momentum eqn at (j mav - 1+1/2) 

„ ITldX 

and for j = 1 

U = 0 

V = 0 (or specified) 

.n-momentum eqn at (1 + 1/2). 

The arrangement of the equations and the boundary 
conditions is shown in Fig. 7. 

For periodic boundaries, the surface normal 
velocities cannot be specified. The total 
number of discrete equations remains 3j max -4. 
However, the number of unknowns is now reduced 
(Fig. 8), since the values of the unknowns at 
j = Jmax are ec l u al to those at j = 1. 

Therefore the total number of unknowns is 
now (3j max - 3). This requires only one addi- 
tional equation to close the system. This condi- 
tion is obtained by applying the 5-momentum 
equation at j = j max . This relies on the fact 
that the location j = j max + 1 is equivalent 
to the location j = 2. The resulting periodic 
block tridiagonal system can once again be 
solved using the standard LU decomposition. 28 

This procedure can be directly extended to 
three-dimensi nal flows since the n and 5 
momentum equations are discretized in exactly 
the same manner as for two-dimensional flows. 

The continuity equation is now dicretized at 
(j - 1/2, k - 1/2) (Fig. 2). Along the lines 
Jmax - 1 and k max - 1 , the arrangement of equa- 
tions i s as fol lows : 

F° r < Jmax - 1 * k) 


‘continuity at <j max - 1 / 2 , k - 1 / 2 )' 

^-momentum at (j mav - 1 , k) 
m& x 

5 -momentum at <j max - 1 , k + 1 / 2 ) 
continuity at <j max - 1 - 1 / 2 , k- 1 / 2 ) 

with a similar form for (j, k max - 1). For inte- 
rior points, 2 < j < j max - 1 , the arrangement 
is given as 

'continuity at (j- 1 / 2 , k- 1 / 2 )' 

5 -momentum at ( j , k) 
n-momentum at (j+l// 2 , k) 

. 5 -momentum at (j, k+ 1 / 2 ) 

The resulting block nine-diagonal matrix 
system is solved in the same manner with the 
MSIP scheme. The boundary conditions can be 
summari zed as fol lows . 

Inflow (5 = 5 0 ): At the inflow boundary, 

velocities are specified. A condition on the 
staggered pressure is required only for the com- 
pressible case (w * 0 ). 

Outflow (5 = 5 max ): the ouflow boundary, 

only one boundary condition on the staggered 
pressure is required. Either p or 9p/35 is 
specified. 

Lower and left wall boundaries (n = 0, and 
5 = 0); No slip and zero injection are specified 
on the solid walls (U = 0, V = 0, and W = 0). 

A boundary condition for pressure is not 
required. The normal momentum equations, at the 
corresponding boundaries, are applied to obtain 
the wall pressure (n-momentum equation at n = 0 
+ An /2 and 5 momentum equation 
at 5 = 0 + A 5 / 2 ) . 

Upper and right boundaries (n = nmax and 
5 = 5max ): Once again zero injection and no slip 

are specified on the solid walls. A boundary 
condition for pressure is not required. The nor- 
mal momentum equations are applied at the bound- 
aries to obtain the surface pressure (n-momemtum 
equation at n = nmax - An/2 and 5-momentum 
equation at 5 - 5max - A5/2). 

Solution Procedure 

Starting from the inflow boundary and then 
at each marching location i, the block nine 
diagonal system shown in Fig. 3 is solved with 
the MSIP algorithm for [U, V, W, P] T . The den- 
sity and temperature for compressible flows, is 
updated after the local iteration for the quasi- 
linearized system has converged at each marching 
location i. The density and temperature are 
evaluated using the state and energy equations. 
The relaxation (marching) procedure proceeds to 
the downstream boundary. The terms with super- 
script (n - 1 ) are then updated from the previ- 
ous marching sweep or global iteration. The 
relaxation process from the inflow to the out- 
flow boundaries is repeated until the maximum 
change in the pressure field for two consecutive 
global iterations is less than a prescribed 
tolerance, e.g. , lO -5 . 
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Stabi 1 i ty 

As discussed in Ref. 19, a detailed stabil- 
ity analysis of the relaxation procedure for two- 
dimensional incompressible flow is presented in 
Ref. 16 and a similar analysis for compressible 
flow is given in Ref. 29. The analysis shows 
that the relaxation procedure for the pressure 
field is unconditionally stable. Since relation 
of the pressure field in the marching direction 
is the same for three- and two-dimensional flows, 
a similar conclusion can be inferred for three- 
dimensional flows. For the cross-plane inver- 
sion, the MSI procedure has been shown to be 
unconditionally stable in Ref. 27 Therefore, 
the overall solution procedure is postulated as 
unconditionally stable. No stability limitations 
were encountered in the present calculations. 

Results 


The three-dimensional flows considered in 
this study were chosen primarily to validate the 
scheme and to compare the results with available 
experimental data and numerical solutions 
obtained by other schemes. First, a simple case 
of laminar developing flow in a straight duct of 
square cross section was considered. As noted 
previously, the Reynolds number is based on the 
uniform inlet velocity and the hydraulic diameter 
H of the cross section. The velocities and 
lengths are nondimensional ized with respect to 
the inlet velocity and the hydraulic diameter 
respectively. The number of grid points used in 
the streamwise direction is 51 and in the cross 
plane 11 by 11. The results obtained for incom- 
pressible flow were compared with the numerical 
solution of Rubin and Khosla, 3 ^ obtained using a 
boundary 1 ayer/potenti al core analysis, and with 
the experimental data of Goldstein and Kreid 31 
(see Figs. 9 and 10). The comparison shows very 
good agreement of the RNS results with both the 
earlier numerical results and the experimental 
data. Next, a slightly more severe case of deve- 
loping flow in a circular arc (curved) duct of 
square cross section was considered (Fig. 11). 

As in the case of the straight duct, the Reynolds 
number is based on the uniform inlet velocity 
and the hydraulic diameter of the cross sec tion . 
The Dean number, which is defined as Re/ \J r 7 h 
is 55, and Reynolds number is 205. Once again, 
only the Incompressible flow solution was 
obtained, as the data for this case was readily 
available. The number of streamwise stations 
was 101, with a grid size of 15 by 15 in the 
cross plane. The development of the streamwise 
velocity profiles in the radial plane and in the 
transverse plane, from the entrance to the fully 
developed region, is shown in Figs. 12 and 13. 

The fully developed streamwise velocity profile 
in the radial plane was compared with the numer- 
ical solutions obtained by Kreskovsky et al., 32 
who assumed a parabolic secondary flow correc- 
tion to the primary flow, Ghia and Sokhey, 9 who 
used a parabolic method and with the experimen- 
tal data of Mori et al. 33 The comparison, as 
seen in Fig. 14, shows good agreement with the 
numerical results. The experimental data dif- 
fers from all of the numerical solutions. This 
disagreement of the experimental results can be 
attributed to possible inaccuracies in the mea- 
surements. Figure 15 shows the comparison of 
the fully developed secondary velocity profile 
in the transverse plane with those predicted by 


Kreskovsky et al. 32 and Ghia and Sokhey. 9 The 
agreement of the RNS solution with the other 
numerical solutions is very good. A vector plot 
of the fully developed secondary velocity in the 
cross section is given in Fig. 16. This clearly 
depicts the plane of symmetry and the location 
of the vorticies that appear away from the axis 
and toward the outer wall due to the effects of 
the centrifugal force. Compressible low Mach 
number solutions were also obtained but are not 
presented here. 


Summary 

A consistent procedure has been developed 
for the application of zero injection and pres- 
sure boundary conditions for the RNS algorithm. 
The matrix structure of the discretized equations 
has been reordered near the boundaries so the 
physically meaningful boundary conditions can be 
applied on the boundaries in a direct manner. 
Global and local mass conservation is satisfied 
automatically without the necessity of iteration 
or a Poisson pressure equation. The modified 
strongly implicit procedure (MSIP) is used to 
invert the nine diagonal matrix resulting from 
the coupled system of equations for velocities 
and pressure in the cross plane. The procedure 
has been validated for two example cases, deve- 
loping flow in straight and curved ducts of 
square cross section. The solutions are in good 
agreement with other computed results and experi- 
mental data. The scheme will be applied to more 
complex geometries of practical interest and to 
compressible flows in future studies. 
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Appendix A - Transformation 


5 = 5<x,y,z) 
n = n(x,y,z) 


C = C(x,y,z) 

The Jacobian J is given by: 


J - 


Wc + Wn + x nV* - Wn “ Wc ' Wc 


The contravariant components U, V, and W, 
written without metric renormalization are given 
by 


U = £ x u + C y v + £ z w 

V = n x u + n y v + n z w 

W = C x u + C y v + C 2 w 

where u, v, and w are the Cartesian velocity 
components in x, y, and z directions respec- 
tively. The shear stress terms appearing in the 
momentum equations are given as follows 
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% 




5 x x xx + ^y x xy + ^z T x 


*z xz 


^x T xy + ^y x yy + ^z x zz 


5 x x xz + ^y x yz + ^z T xz 


" Vxx + n y x yx + n z T xz 


n x T xy + Vyy + n z x yz 


T z = n x x xz + n y x yz + n z x zz 


■ Vxx + S T y* + ^ zTzx 


^x x xy + ^y x yy + ^z x zy 


T z = ^x x xz + ^y x yz + ^z x zz 

where t xx , x X y, etc. are the regular Cartesian 
shear stress components. The Cartesian deriva- 
tives are expanded in n, and £ space via 
chain rule relations such as u x = 5 x u^ + nx u n 
+ £ x ur. The Cartesian velocity components are 
in turn expressed in terms of the contravari ant 
velocity components defined earlier by the 
fol lowi ng 

u = x r U + x V + x r W 
£ n 

v = y^U + y n V + y ? W 

w = z-U + z V + z r W 
% n <> 

Finally, the terms G^, q^ y etc. appearing in 
the K are given as 

G CC = U(<5 xV + (5 yV + ( Vs W> 

+ V((5 x ) n u + (5 y ) n v + (5 z ) n w) 

+ M((5 K ) t u + <5 y ) ; v + <? z ) ; w) 

- U((n x )^u + (n y )^v + (n 2 )|W) 

+ V((n x ) n u + <n y ) n v = <n z > n w> 

+ W((n x >^u + (n y )^v + <n 2 >^w) 

G CC = U(< Vc u + ( Vc v + <C zV > 

+ V(< V n u + ( Vn v + <t z > n w) 

+ H((C x > c u + <C y > c v + (C 2 > c w) 

= 5 X 2 ♦ 5 y 2 ♦ 5 Z 2 


ORIGINAL PAGE IS 
OF POOR QUALITY, 
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,ne 


= C X C X ♦ C y S y ♦ S 2 < z 


5 x n x + ? y n y + C 2 n 2 


g 

g« 

9 Cn 


nn 2 2 2 

‘ n x + n y + n z 
= n x C x ♦ n y S y + n z C 2 


“ <x«x + W ^ 


“ ^x n x + S n y + c z n z 


= Vx + SS + ? z<z 
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FIGURE 1. - DISCRETIZATION IN THE ^-DIRECTION. 
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FIGURE 2. - DISCRETIZATION IN THE CROSS-PLANE. 
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FIGURE 4. - COMPUTATIONAL MOLECULE IN THE CROSS PLANE. 



FIGURE 5. - DISCRETIZATION FOR 
2-D CASE IN THE NORMAL 01) 
DIRECTION. 
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FIGURE 6. - BLOCK TRIDIAGONAL SYSTEM AT MARCHING LOCATION (i) FOR 2-D FLOW. 



FIGURE 7. - GROUPING OF DIS- 
CRETIZATION GOVERNING EQUA- 
TIONS IN THE FLOW FIELD. 
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FIGURE 8. - PERIODIC BOUNDARIES 
FOR 2-D CASE. 
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CENTERLINE VELOCITY 



FIGURE 9. - CENTERLINE VELOCITY VARIATION. 



FIGURE 10. - AXIAL VELOCITY PROFILE AT 
£ = 0.295, Z = 0. 
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FIGURE 11. - CURVED DUCT GEOMETRY AND COORDINATE SYSTEM. 



FIGURE 12. - DEVELOPMENT OF PRIMARY VELOCITY PROFILE FOR 
A CIRCULAR ARC DUCT OF SQUARE CROSS SECTION, Re = 205, 
R/H = 14, K = 55. 
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FIGURE 13. - DEVELOPMENT OF PRIMARY VELOCITY PROFILE 
IN THE TRANSVERSE PLANE FOR A CIRCULAR ARC DUCT OF 
SQUARE CROSS SECTION, Re = 205, R/H = 14, K = 55. 
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FIGURE 14. - FULLY DEVELOPED PRIMARY FLOW VELOCITY 
PROFILE FOR A CIRCULAR ARC DUCT OF SQUARE CROSS 
SECTION, R/H = 14, Re - 205, K = 55. 
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FIGURE 15. - FULLY DEVELOPED RADIAL VELOCITY 
PROFILE AT y/H = 0.4 FOR A CIRCULAR ARC 
DUCT OF SQUARE CROSS SECTION, R/H = 14, 

Re = 205, K = 55. 
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CIRCULAR ARC DUCT. Re = 205, R/H = 14, K = 55. 
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